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A  ONE-DIMENSIONAL  TIME -DEPENDENT  MODEL 
FOR  FLAME  INITIATION,  PROPAGATION  AND  QUENCHING 

1.  INTRODUCTION 

This  report  describes  a  one-dimensional,  time -dependent,  Lagrangian 
numerical  model  developed  to  study  the  initiation,  propagation  and  quenching 
of  laminar  flames.  The  model  incorporates  a  number  of  new  approaches  and 
algorithms  which  have  now  been  tested  by  conparisons  to  less  complex  or 
analytic  solutions  and  by  comparisons  to  experimental  data.  These  new  ele¬ 
ments  include:  ADINC  (l)  an  inplicit,  lagrangian  method  for  solving  the 
convective  parts  of  the  conservation  equations;  DFLUX  (2,3)  a  variable 
accuracy  algorithm  for  determining  diffusion  fluxes  without  having  to  invert 
matrices;  VSAIM,  a  vectorized  version  of  the  ordinary  differential  equation 
solver,  CHEMEQ  [U,5l;  and  a  new  method  for  treating  an  open  boundary  in  an 
implicit,  Lagrangian  calculation.  An  asymptotic  coupling  method  is  used  in 
conjunction  with  timestep  splitting  to  couple  the  various  physical  and  chemi¬ 
cal  processes.  ‘Ibis  approach  allows  the  use  of  entirely  different  algorithms 
for  the  physical  processes  represented  by  the  different  terms  in  the  conser¬ 
vation  equations. 

Manuscript  submitted  July  26, 1982. 


1 


The  model  has  "been  used  for  a  variety  of  flame  studies  of  hydrogen- 
oxygen-nitrogen  mixtures.  These  include  calculations  of  minimum  ignition 
energies  (6,71,  flammability  limits  [8],  quench  volumes  (7l»  and  burning 
velocities  (9].  The  chemical  rate  scheme  has  now  been  tested  extensively,  as 
have  the  input  parameters  for  the  transport  properties.  Thus  we  expect  the 
model  to  calculate  correctly  the  time-dependent  behavior  inplied  by  the 
initial  and  boundary  conditions  supplied. 

After  a  brief  review  of  related  literature,  a  summary  is  given  of  the 
inportant  numerical  features  of  the  model.  Then  several  calculations  are 
presented.  The  first  is  a  flame  initiation  and  minimum  ignition  energy  study 
of  a  mixture  of  H2:02:N2/2:1:10  initially  at  298  K  and  1  atm.  This  is 
inherently  a  time-dependent  problem  and  takes  advantage  of  this  property  of 
the  model.  The  second  problem  discussed  is  a  calculation  of  the  burning 
velocity  of  the  mixture  H2:02:N2/2:1:4,  again  at  298  K  and  1  atm.  The 
burning  velocity  describes  a  steady  state  property  of  the  system;  thu3  we 
show  below  that  the  model  also  does  well  on  this  type  of  calculation. 
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2.  BACKGROUND 


The  laminar  flame  problem  is  the  earliest  combustion  problem  to  be 
studied  theoretically  which  required  the  simultaneous  consideration  of  both 
f  fluid  dynamics  and  chemical  kinetics  for  its  solution  llO].  The  problem  of 
determining  the  propagation  velocity  of  a  deflagation  wave  was  first  studied 
by  Mallard  and  LeChatelier  (ill.  Since  then  there  have  been  many  analytical 
attempts,  [e.g.  12,  13l  with  varying  degrees  of  complexity  and  approximation, 
to  study  simple  steady  flames.  Some  of  these  have  been  summarized  by 
Williams  [lOl .  Whereas  these  analytical  techniques  are  extremely  informa¬ 
tive,  they  cannot  satisfactorily  describe  the  detailed  structure  of  flames 
since  they  all  include  very  approximate  chemical  kinetic  schemes  and  trans¬ 
port  properties.  For  this  one  has  to  resort  to  numerical  methods. 

One  of  the  first  attempts  to  solve  the  laminar  flame  problem  numerically 
was  that  of  Hirschfelder  and  co-workers  [l4].  They  formulated  the  unsteady 
flame  problem  as  a  system  of  three-dimensional  nonlinear  partial  differential 
equations  and  they  solved  the  one-dimensional  steady  flame  as  a  two  point 
boundary  value  problem.  They  used  approximate  methods  to  estimate  the  mass 
flow  rate  which  is  the  eigenvalue  of  the  two-point  boundary  value  problem  and 
then  used  a  numerical  shooting  method  to  obtain  the  temperature  and  species 
,  profiles.  Although  the  first  problem  they  studied  involved  only  3ingle  step 
kinetics,  they  later  applied  the  same  solution  procedure  to  study  flames 

V 

for  which  the  kinetics  involved  chain  reactions  (l5»l6J. 

In  1956,  Spalding  (lTl  applied  a  time-dependent  system  of  nonlinear 
partial  differential  equations  using  explicit  finite-difference  techniques  to 
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study  hydrazine  flames.  He  assumed  initial  profiles  for  the  temperature  and 
species  concentrations  and  obtained  the  steady  state  burning  velocity  by 
carrying  out  the  confutations  for  a  sufficiently  long  time.  Spalding's 
method  was  subsequently  used  by  Adams  and  Cook  [18]  to  study  the  effect  of 
pressure  on  the  reaction  mechanism  and  speed  of  hydrazine  flames,  and  by 
Dixon-Levis  (19]  to  study  rich  hydrogen-oxygen  flames.  DLxon-Levis  [19)  used 
a  set  of  fourteen  reactions  to  describe  hydrogen-oiygen  kinetics,  but  assumed 
a  steady  state  approximation  for  the  radical  distribution  in  the  flame.  In 
1968,  Dixon-Levis  (20)  studied  the  flame  structure  and  flame  reaction  kinet¬ 
ics  of  hydrogen-oxygen-nitrogen  flames  using  detailed  expressions  for  the 
diffusive  transport  coefficients  in  railticomponent  systems.  He  also  used 
this  time -dependent  numerical  model  to  stu<ty  ignition  by  localized  sources 

[ail  • 

Hie  work  of  Dixon-Levis  and  co-workers  described  above  proved  the  one¬ 
dimensional  unsteady  flame  with  coniplex  chemistry  and  detailed  diffusive 
transport  coefficients  could  be  solved  using  numerical  methods.  Hie  emphasis 
then  shifted  toward  the  development  of  more  economical  and  iuproved  numerical 
methods.  In  1971,  Spalding,  et  al*  [22]  presented  an  implicit-finite  differ¬ 
ence  method  in  which  the  unsteady  flame  equations  were  transformed  to  a  form 
which  could  be  solved  by  a  numerical  method  developed  by  Patankar  and 
Spalding  [23]  for  two  dimensional  boundary  layer  equations.  They  used  sim¬ 
plified  transport  properties  [14]  and  a  four  step  chain  reaction  mechanism  to 
study  the  propagation  of  hydrogen-bromine  flames  [24]..  However  this  method 
requires  careful  ordering  of  the  solution  of  the  equations,  and  a  certain 
manner  of  linearizing  the  source  terms.  Hie  same  numerical  procedure  was 
used  later  by  Tsatsaronis  [25]  to  study  unsteady  flame  propagation  in 
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methane,  oxygen,  nitrogen  mixtures  using  very  detailed  chemical  kinetic 
mechanisms  and  transport  properties* 


The  kinetics  and  propagation  of  methane-air  flames  was  studied  by  Smoot 
et  al.  (26]  using  a  mixed  explicit-implicit  finite  difference  technique.  Hie 
diffusive  transport  terms  were  solved  explicitly  and  the  Kinetic  terms  were 
solved  using  linearized  implicit  techniques.  Bledgian  [27]  employed  a  method 
of  lines  technique  where  the  nonlinear  parabolic  initial-boundary  value 
problem  is  reduced  to  a  set  of  nonlinear  first  order  initial  value  problems. 
Sophisticated  initial  value  problem  integrators  were  used  to  solve  the 
resulting  ordinary  differential  equations  system.  Instead  of  using  finite 
difference  approximations  to  discretize  the  spatial  derivative,  Margolis  [28) 
used  a  spline  collocation  procedure  in  his  method  of  lines  approach. 

Westbrook  and  Dryer  (29)  have  deduced  a  comprehensive  reaction  mechanism  for 
methanol  oxidation  involving  eighty  four  reactions  using  an  illicit  finite- 
difference  algorithm  (30)  to  solve  the  unsteady  flame  equations.  They  used 
sinplified  expressions  for  the  diffusive  transport  coefficients  which  were 
adjusted  to  give  good  laminar  flame  speed  prediction  for  methane-air  flames 

Dixon-Lewis  has  more  recently  used  a  conposite  flux  method  [3l]  to  study 
a  variety  of  problems  like  the  kinetic  mechanism,  structure  and  propagation 
of  flames  in  hydrogen-oxygen-nitrogen  flames  [32]  and  flame  inhibition  by 
organic  halogen  compounds  (33l •  He  discusses  the  ranges  of  applicability  of 
the  partial  equilibrium  and  quasi-steady  assunptions  in  relation  to  the  dis¬ 
tribution  of  radical  populations  in  the  flames  [ 32) . 

Warnatz  has  extensively  studied  both  freely  propagating  and  burner 
stabilized  flames  in  a  variety  of  premixed  gases  [3^,  35>  36,  37,  38).  He 
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linearizes  the  chemical  reaction  terms  and  solves  the  time-dependent  equa¬ 
tions  illicitly  with  assumed  Initial  guesses  for  the  temperature  and  species 
profiles.  He  also  uses  a  simplified  transport  model  [3M  which  agrees  well 
with  the  complete  nulticonponent  fornulation  of  diffusion  and  thermal  conduc¬ 
tion.  In  his  recent  study  (38]  of  the  concentration,  pressure  and  tempera¬ 
ture  dependence  of  the  flame  velocity  in  hydrogen-oxygen-nitrogen  mixtures  he 
concludes  that  at  the  present  state  of  knowledge  predictions  of  laminar, 
premixed  flame  propagation  should  be  as  reliable  as  experimental  results. 

Coffee  and  Heimerl,  using  a  method  of  lines  approach  [39]  examined  dif¬ 
ferent  methods  of  approximating  niiltispecies  transport  phenomena  in  models  of 
premixed,  laminar,  steady -state  flames.  They  concluded  that  the  selection  of 
tne  input  values  for  the  individual  species  transport  properties  is  more 
important  than  the  selection  of  the  approximation  method  [40] . 

There  have  also  been  a  number  of  approaches  which  solve  a  steady  state 
fornulation  of  the  flame  equations.  These  include  the  works  of  Dixon-Lewis 
(Ull ,  Wilde  [42]  and  Sraooke  [43] •  The  advantages  and  the  difficulties  in 
using  steady  state  solution  procedures  have  been  discussed  by  Smooke  [43] . 
These  methods  are  good  for  obtaining  burning  velocities  and  steady  state 
profiles,  but  cannot  provide  information  on  the  ignition  and  development  of 
flames. 

The  model  presented  in  this  report  solves  the  time-dependent  partial 
differential  equations  and  hence  can  describe  the  initiation,  propagation  and 
quenching  of  laminar  flames.  It  uses  detailed  chemical  kinetic  models 
without  any  steady  state  or  quasi -equilibrium  assunptions  and  so  can  be  used 
to  study  all  laminar  premixed  flames.  It  uses  an  asymptotic  coupling  method 
in  conjunction  with  timestep  splitting  to  couple  the  various  physical  and 
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chemical  processes,  This  approach  allows  the  use  of  entirely  different  algo¬ 
rithms  for  the  processes  represented  by  the  different  terms  in  the  conserva¬ 
tion  equations.  Thus  it  is  numerically  very  efficient  and  inexpensive.  A 
Lagrangian  formilation  is  used  in  order  to  maintain  steep  gradients  for  a 
long  period  of  time.  The  model  has  been  used  for  a  variety  of  flame  studies 
of  hydrogen-oxygen-nitrogen  mixtures.  'Uiese  include  calculations  of  ndninum 
ignition  energy  (6,71 »  flammability  limits  [81,  quench  volumes  (?1»  and 
burning  veloc it i es  ( 9 1 • 
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3.  NUMERICAL  MODEL 


3.1  Basic  Equations 

We  solve  the  time-dependent  equations  for  conservation  of  total  mass 
density  p,  momentum  pv,  and  energy  density  E^  as  well  as  the  individual 
species  number  densities  {nj}.  These  may  be  written  as  [2,10] : 

Jt  *  p-  (1) 


It  ■  *  pj  -  Vi 


•gY'q  -7.(  pw)  -  VP  +  V.  7v  +  (Vv)T| 


*  -V.Ev  -  V.Pv  -  V.Q 


where  the  heat  flux  vector,  Q,  is  defined  as 


Q  =  -XV  T  +  In.h.V.  +  F l  kJ  V 
”  a-  |  3  J1  j  J  - 


The  quantity  v  is  the  fluid  velocity,  the  (Vj }  are  the  diffusion  veloci¬ 
ties,  and  {Pj  }  and  {Lj }  refer  to  the  chemical  production  and  loss  pro¬ 
cesses  for  the  individual  species  j.  The  quantities  nm  and  Am  are  the 
mixture  viscosities  and  thermal  conductivities  respectively.  The  superscript 
"T"  in  the  last  term  of  Eq.  (3)  indicates  that  the  transpose  is  taken.  The 
quantity  E  in  Eq.  {k)  is  the  total  energy  Ep  minus  the  total  heat  of  forma¬ 
tion.  The  quantities  {hj }  are  the  temperature  dependent  enthalpies  for 
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each  species,  and  the  {Djk}  and  {KTj }  are  Vie  sets  of  binary  diffu¬ 
sion  coefficients  and  thermal  diffusion  ratios,  respectively. 

We  also  assume  that  the  mixture  behaves  like  an  ideal  gas  so  that  the 
pressure,  P,  may  be  written  as: 


P  »  NkgT  (6) 

where  N  is  the  total  number  density,  kg  is  Boltzmann’s  constant  and  T  is 
the  temperature.  The  model  presented  in  this  report,  however,  is  not  re¬ 
stricted  to  ideal  gases  and  in  fact  any  equation  of  state  may  be  used. 

The  diffusion  equations  may  be  written  as 


(7) 


where  the  source  terms  Sj 


S 


J 


are  defined  as 


(8) 


The  diffusion  velocities  are  also  subject  to  the  constraint: 


I ' 0 


(9) 


Below  we  describe  the  basic  features  of  the  solution  techniques  adopted 


to  solve  the  terms  representing  the  chemical  reactions,  the  convective  and 
diffusive  transport  processes,  and  techniques  for  coupling  the  solutions  of 
the  equations  representing  these  processes. 


3.2  Convective  Transport 

The  convective  transport  terns  in  Eqs.  (1-4)  are  solved  by  the  algorithm 
ADINC  [l]  .  ADINC,  designed  for  either  Adiabatic  or  incompressible  flows,  is 
an  illicit,  Lagrangian  algorithm.  Since  ADINC  communicates  compression  and 
expansion  across  the  system  illicitly,  it  overcomes  the  Courant  time-step 
limit.  Since  it  is  Lagrangian,  it  can  maintain  steep  gradients  computation¬ 
ally  for  a  long  period  of  time.  This  is  important  in  flame  calculations 
where  the  diffusive  transport  of  material  and  energy  can  govern  the  system 
evolution  and  therefore  must  be  calculated  accurately. 

ADINC  solves  the  following  equations  for  mass  and  momentum  transport  in 
one  dimension: 

fg-  *  -  P  V.v  (10) 

dv 

p  55  -  F  (n) 

The  energy  evolution  equation  is  eliminated  by  using  an  adiabatic  equation  of 
state: 

p(P,S)  *  Pc  +  (P/S)1/Y  (12) 

This  equation  of  state  with  pc  *  0  is  correct  for  adiabatic  compression  and 
expansion  of  an  ideal  gas.  The  entropy,  S,  is  assumed  constant  throughout 
the  numerical  integration.  Non-adia'oatic  processes,  such  as  external 
heating,  thermal  conduction  and  chemical  energy  release,  are  added  to  Eqs. 

( 10-11)  using  the  timestep  splitting  methods  described  below. 

Given  an  aproximation  to  the  pressure,  ADINC  calculates  the  fluid  den¬ 
sity  using  Eq.  (12).  This  equation-of-state  density  is  compared  to  the  den¬ 
sity  derived  from  the  fluid  dynamics  through  Eq.  (10).  The  difference  is 


10 


iterated  to  zero  using  a  quadrat ically  convergent  implicit  solution  of  Eq. 

(ll)  which  then  gives  an  inproved  approximation  to  the  pressure.  During  this 
iteration  the  analytic  derivative  dA/dP  is  used  where  A  is  the  volume  of  a 
computational  cell.  Thus 

xIf-  (P/s)1/T  <«> 

for  the  particular  equation  of  state  (12). 

ADINC  assumes  that  pressure  and  density  are  constant  within  each  indi¬ 
vidual  finite-difference  cell  and  that  the  physics  is  evolving  slowly  enough 
for  full  communication  across  that  cell  to  have  occurred  in  a  timestep.  When 
the  densities  in  adjacent  cells  are  different,  the  pressure  gradient  between 
the  two  cells  will  impart  a  different  acceleration  to  the  fluid  Just  to  the 
right  and  to  the  left  of  the  interface  between  the  two  cells.  If  the  fluid 
were  permitted  to  move  according  to  these  distinct  accelerations,  the  fluid 
would  either  overlap  or  a  gap  would  appear  at  the  interface  after  a  short 
while.  To  prevent  this,  ADINC  matches  the  acceleration  calculated  from  the 
left  to  that  from  the  right.  This  is  described  in  more  detail  by  Boris  [l]. 

Problems  can  be  set  up  in  one  of  four  geometries:  cartesian  coordi¬ 
nates,  cylindrical  coordinates,  cylindrical  coordinates,  spherical  coordi¬ 
nates  and  a  variable  power  series  coordinates  for  treating  one-dimensional 
nozzle  like  geometries.  The  boundary  conditions  are  controlled  by  specifying 
the  position  and  velocity  of  the  first  and  last  cell  interface. 


3«3  Chemical  Kinetics  Calculations 


The  coupled,  nonlinear,  ordinary  differential  equations  vhich  describe 
the  chemical  interactions  are  taken  from  that  part  of  Eq.  (2)  which  repre¬ 


sents  the  production  and  loss  of  reacting  species: 

-  PJ  -  LJ  nj  • J 


3t 


.  M, 


(1U) 


where  M  is  the  total  number  of  species  present.  The  functional  dependencies 
of  the  terms 

pj  -  PjK(t)) 


Lj  «  Lj{nk(t)}  ,  k  -  1,  M  (15) 

emphasize  the  strong  coupling  between  the  various  species.  The  system  repre¬ 
sented  by  Eq.  (l4)  may  be  stiff  when  there  are  large  differences  in  the  time 
constants  associated  with  different  chemical  reations.  Stiffness  may  occur 
for  different  species,  in  different  locations,  at  different  times  or  simul¬ 
taneously  throughout  the  course  of  an  integration.  Because  of  this,  when 
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there  are  a  large  number  of  reactions,  the  solution  of  the  chemical  kinetics 
equations  is  usually  the  most  expensive  part  of  a  detailed  reactive  flow 
calculation.  Furthermore,  the  computational  cost  increases  with  the  number 
of  species  and  the  dimensionality  of  the  problem.  Therefore  a  method  which 
is  efficient,  accurate,  conservative,  stable  and  which  does  not  require  stor¬ 
age  of  large  quantities  of  data  from  one  timestep  to  another  is  required. 

Such  a  method  is  VSAIM,  which  is  a  fully  vectorized  version  of  the  selected 
asymptotic  integration  method  employed  in  CHEMEQ  (4,5). 


In  VSAIM,  the  stiff  equations  are  identified  and  solved  using  a  very 
stable  asyuptotic  method  while  the  remaining  equations  are  solved  using  a 
standard  classical  method*  Since  it  is  self  starting  and  does  not  require 
matrix  inversion,  it  is  inexpensive.  Species  densities  at  only  one  time 
level  are  needed.  Accuracy  may  be  controlled  by  predetermined  convergence 
parameters  which  control  the  stiffness  criterion,  the  timestep,  and  thus 
effectively  the  degree  of  conservation.  Hie  method  is  stable  for  large  time- 
steps  and  is  vectorized  for  use  in  parallel  processing  computers. 


3.4  Diffusive  Transport 

The  diffusive  transport  processes  considered  in  this  model  are  molecular 
diffusion,  thermal  diffusion  and  thermal  conduction.  These  are  the  parts  of 


Eqs.  (2,4)  which  are  represented  by  the  expressions: 

an, 


at 


-  v.[-xvt  +  Xh^x,  +  pIk^  v1 ) 


J  i-i 


J  -J 


(16) 

(17) 


J  °  °  u  1 

These  processes  are  crucial  to  the  description  of  flame  ignition  and  propaga¬ 
tion  since  they  are  the  mechanism  by  which  heat  and  reactive  species  are 
transported  ahead  to  the  unburned  gas.  The  effects  of  viscosity  have  not 
been  considered  in  this  paper. 

The  diffusion  velocities  {Vj }  in  Eqs.  (16,  17)  are  determined  by 
solving  the  equations  (2,44] 


M 

=  l 


!4Mv-v) 


7  (nj/N)  +  Kj  f 


(18) 
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(19) 


subject  to  the  constraints 

M  M 

l  S  »  0,  and  l  p  V  «  0. 

J«1  J  j-1  0  J 

An  iterative  dgorithra,  DFLUX,  enables  us  to  obtain  the  diffusion  velocities 

without  the  cost  of  performing  matrix  inversions  (2,3 ] •  Ibis  method  is  of 

0(M2)  and  is  vectorized.  Ihus  it  is  substantially  faster  than  0(M3)  matrix 

inversions  when  four  or  more  species  are  involved. 

The  Eqs.  (l6,17)  are  conservatively  differenced  and  solvsd  explicitly. 

Solving  these  equations  requires  knowledge  of  the  binary  diffusion  coeffi- 

T 

cients  Dj^,  the  thermal  diffusion  ratios  Kj  and  the  thermal 
conductivity  X.  Representations  of  these  coefficients  which  are  easily  used 
in  reactive  flow  models  and  do  not  require  the  expensive  inversion  of  matri¬ 
ces  are  chosen  and  these  are  discussed  in  a  subsequent  section  of  this 
report. 

3» 5  Timestep  Splitting 

In  the  asymptotic  timestep  split  approach  (2) ,  the  individual  terms  in 
the  Eqs.  (1-4)  are  solved  independently  as  described  above  and  then  aaytqptot- 
ically  coupled  together.  Since  both  chemistry  and  sound  waves  are  usually 
stiff  in  deflagration  problems,  special  care  is  required  in  coupling  the 
chemical  heat  release  to  the  fluid  dynamics.  In  a  flame,  fluid  dynamic 
expansion  and  diffusive  transport  relieve  the  pressure  from  the  flame  region 
as  fast  as  it  is  generated.  Thus  the  pressure  stays  effectively  constant. 
Small  pressure  fluctuations,  0(v.2/c2),  do  exist  and  are  Just  large  enough  to 
drive  the  flows  which  reapportion  the  energy  released  by  chemistry  or  trans¬ 
ported  by  diffusive  processes. 
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The  energy  released  in  chemcal  reactions  and  that  transported  by  diffu¬ 
sion  and  thermal  conduction  are  converted  to  effective  pressure  changes  as 
described  below.  Before  the  fluid  dynamic  timestep,  the  sum  of  all  these 
pressure  changes  is  added  to  the  old  pressure  distribution  to  obtain  a  new 
pressure  distribution.  This  new  pressure  distribution  is  used  to  modify  the 
cell  entropies.  The  new  pressure  distribution  and  cell  entropies  are  used  in 
the  convection  transport  algorithm,  ADINC  to  obtain  the  convect  fluid  expan¬ 
sion  as  described  earlier. 

The  chemistry  step  should  be  taken  at  constant  pressure,  but  if  the 
profiles  change  slowly  per  timestep,  it  may  also  be  taken  at  constant  volume 
with  temperature  held  fixed.  On  completion  of  the  chemistry  integration,  the 
heat  release  is  converted  to  an  effective  pressure  change  at  constant  volume 
because  the  cell  volume  has  been  held  fixed  during  the  chemical  kinetics 
calculation.  The  chemical  energy  released  can  be  calculated  from  the  heats 
of  formation  and  the  changes  in  the  number  densities  of  the  various  species. 
Defining  an  energy  density  e  by: 

*  -  et  -  -  \  Vj  (20) 

J 


at  the  end  of  each  chemical  time  step  we  have 

e(t+At)  *  e(t)  +  I  (n.(t+At)  -  n,(t))h  . 

JJ  0  oj 


The  changes  in  the  number  densities  are  calculated  by  solving  the  chemical 


rate  equations  as  described  in  an  earlier  section  of  this  report.  Using  the 


definition  of  enthalpy,  the  ener®-  density  e  is  written  as: 


c  =  l  h j(T)rij  -  P. 

J 


The  temperature  at  the  end  of  the  chemical  time  step  is  evaluated  from 


Eq.  (22)  using  the  new  values  for  e  and  nj  and  the  temperature  dependence 
of  the  sensible  enthalpy  hj.  Since  Eq.  (22)  is  highly  nonlinear  in  temper¬ 
ature,  a  Nevton-Raphson  iteration  technique  is  used  and  it  usually  converges 
in  one  or  two  iterations.  Hiis  is  the  temperature  that  would  have  been 
attained  if  the  entire  energy  released  goes  into  heating  the  fluid.  Since 
the  chemical  energy  released  also  goes  into  fluid  expansion,  it  would  be 
incorrect  to  use  this  as  the  temperature  for  the  next  timestep.  Therefore, 
the  energy  released  is  converted  to  an  effective  pressure  change  as 

1P  ■  V  }  “j  -  Pold  (23) 

where  T  is  the  temperature  calculated  from  Eq.  (22)  and  n^  are  the  new 
number  densities  of  the  various  species.  The  energy  transported  by  thermal 
conduction  and  diffusion  are  also  converted  to  effective  pressure  changes  by 
using  the  equation 

AP  *  AE  (y-1)  (24) 

These  pressure  changes  are  used  as  energy  source  terms  for  the  fluid  dynamic 
timestep  as  described  above. 

3.6  Open  Boundary  Condition 

In  the  study  of  unconfined  flames,  a  representation  for  the  open¬ 
boundary  is  required  since  the  size  of  the  conputational  domain  is  limited. 
One  approach  is  to  allow  the  computational  cells  to  increase  in  size  as  they 
get  further  from  the  flame.  Thus  the  computational  domain  is  very  large  and 
there  is  no  corresponding  increase  in  computer  storage.  However,  the  cell 
stretching  should  be  gradual  to  limit  inaccuracies  which  arise  as  a  result  of 
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the  varying  cell  size.  Furthermore,  adaptive  gridding  will  be  essential  with 
this  approach,  since  as  the  flame  front  moves  toward  the  large  cells,  the 
resolution  will  have  to  decrease  around  the  flame. 

A  second  approach  which  ve  have  used  to  study  the  propagation  of  uncon¬ 
fined  flames  takes  advantage  of  the  observation  that  the  effect  of  an  open 
boundary  is  to  keep  the  pressure  essentially  constant.  Therefore,  in  this 
approach  the  movement  of  the  cell  boundary  that  is  on  the  open  end  is  calcu¬ 
lated  by  allowing  the  pressure  in  each  cell  at  the  end  of  a  timestep  to  relax 
adiabatically  to  the  pressure  before  the  timestep.  In  practice  this  is  done 

by  calculating  the  changes  in  the  cell  volume  V  ,  hy 

fP  11/t 

Vc(t+«)  -  Vc(t)[pffyJ  (25) 

where 

«  P(t)  +  AP  (26) 

new 

and  (AP)  is  the  change  in  the  presure  due  to  diffusive  transport,  energy 
released  in  chemical  reactions  and  any  energy  addition  to  the  system.  The 
changes  in  the  volume  of  the  cells  causes  the  location  of  the  open  boundary 
to  change.  The  location  of  the  open  boundary  and  the  fluid  velocity  in  the 
last  cell  (which  is  also  the  velocity  of  the  open  boundary)  are  used  as  open 
boundary  conditions  to  the  convective  transport  algorithm  ADINC. 
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4.  EVALUATION  OF  MATERIAL  PROPERTIES 
In  order  to  apply  the  numerical  model  described  above,  ve  need  expres¬ 
sions  for  the  diffusive  transport  coefficients,  the  enthalpies  of  the  various 
species  involved,  the  chemical  kinetic  rate  scheme  and  the  reaction  rate 
constants.  In  this  section  ve  first  present  general  equations  for  the  evalu¬ 
ation  of  some  of  these  material  properties  and  then  present  the  input  param¬ 
eters  used  for  the  study  of  pre-aixed  hydrogen-oxygen-nitrogen  systems. 


4.1  General  Expressions 

The  diffusive  transport  coefficients  presented  here  have  been  taken  with 
some  modification  from  the  summary  by  Picone  (U5 1 .  The  forms  given  are  for 


mixtures  of  neutral  gases.  Their  derivations  are  based  on  or  are  extensions 
of  the  fundamental  work  of  Chapman  and  Covling  [U4]  and  Hirschfelder,  Curtis 
and  Bird  [46] . 


4.1.1  Thermal  Conductivity 


For  the  mixture  thermal  conductivity,  an  extremely  useful  equation  has 


been  given  by  Mason  and  Saxena  [47]: 

Xm  *  +  2^2  nj  (27) 

where  Xj  is  the  thermal  conductivity  of  a  pure  polyatomic  gas  of  species  j, 
and  is  given  by 

r  X?  M,  l. 


T  *1  M4  1: 

[1  +  (^)l/2  (^)1A 

M.ll/2 

L 
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where  X9  is  the  thermal  conductivity  of  gas  J  if  it  is  monatomic,  or 
J 

it  is  the  thermal  conductivity  of  a  polyatomic  gas  with  internal  degrees  of 
freedom  "frozen'',  and  Mj  is  the  atomic  mass  of  species  J.  She  thermal 
conductivity,  Xj  of  a  pure  polyatomic  gas  can  be  estimated  from  the  mon¬ 
atomic  thermal  conductivities  by  using  the  modified  Euken  factor 
(48),  Ej»  a3 

w;  <29> 


where 


Ej  «  0.115  +  0.354 

B 


(30) 


and  Cpj  is  the  molecular  specific  heat  at  constant  pressure.  The  {X0^} 
may  be  written  as 

L/2 

•1 
‘i 


,n  _  8.322  x  103  rT  ,1' 

XJ  (275)-*“  (m7' 

iii 


(31) 


2)* 

■-fhere  '  ,  a  collision  integral  normalized  to  its  rigid  sphere  value  is 

£ 

a  function  of  T^,  is  the  reduced  temperature  given  by 


(32) 


The  terms  Oj  and  are  force  constants  in  the  potential  function  which 
describes  the  interaction  between  two  molecules  of  the  same  species,  j. 


4.1.2  Ordinairy  Diffusion  Coefficients 

We  assume  that  the  "ordinary"  (concentration)  diffusion  coefficients  for 
a  pair  of  species  (j,k)  in  a  mlticonponent  mixture  is  to  a  first  approxima¬ 
tion  equal  to  that  in  a  binary  mixture  of  species  J  and  k.  For  some  binary 
mixtures  of  dilute  gases,  Mason  and  Marrero  [491  have  given  a  semi-enpirical 
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expression  which  describes  the  variation  of  Dj^  over  a  wide  temperature 
range.  Their  expression  for  Djk  can  be  written  as 

A  T  ^ 

When  the  coefficients  Aj^  and  are  not  available,  the  may 


(33) 


be  estimated  from 


n  2.628  x  10-3  (T3(Mj  + 
Jk  ** 


m 


(1 1)* 

where  ,  a  collision  integral  normalized  to  its  rigid  sphere  value  is 


a  function  of  the  reduced  temperature 


T* 

Jk 


V 


(35) 


The  factors  Oj^  and  Cjk  in  the  above  equation  are  force  constants  in 
the  potential  function  which  describes  the  interactions  between  molecules  of 
species  j  and  k.  In  most  cases  values  of  and  are  not  avail¬ 
able  and  are  usually  calculated  from 


and 


V  -  2  (oj  +  V 


EJt:  -  tejek)1/2' 


(36) 


(37) 


4.1.3  Thermal  Diffusion  Ratios 

Thermal  diffusion  is  a  second  order  effect  (44,  46]  which  is  only  of 
importance  when  there  are  large  differences  between  the  atomic  masses  of 
constituent  species.  An  expression  for  the  thermal  diffusion  ratio,  which  is 


a  measure  of  the  relative  importance  of  the  thermal  to  the  ordinary  diffu¬ 
sion,  can  be  written  as  [44] 


K 


5^ 


(6ck  - »  ,Wk  -  Wj, 


J* 


M,  + 


"k 


(38) 


where  aj  is  the  contribution  of  the  Jth  species  to  the  mixture  thermal 
conductivity  (assuming  that  the  internal  degrees  of  freedom  are  frozen): 


a 


J 


‘j !l  + 


1.065 


i 

k*J 


nk*jk 


-1 


(39) 


and  C*^  is  a  ratio  of  collision  integrals 


(40) 


rn 

The  expression  for  K  j  has  been  derived  by  assuming  that  the  gases 
in  the  mixture  are  monatomic  or  that  the  internal  degrees  of  freedom  are 
frozen.  However,  Chapman  and  Cowling  [44]  have  conjectured  that  the  internal 
degrees  of  freedom  of  the  gas  molecules  is  expected  to  have  a  smaller  effect 
on  the  thermal  diffusion  ratio  than  on  the  thermal  conductivity. 


4.1.4  Chemical  Kinetics 

The  chemical  kinetic  rate  scheme  is  very  specific  for  the  problem  being 
solved  and  the  scheme  [ 50]  used  for  the  hydrogen-oxygen-nitrogen  system  is 
presented  later.  The  reaction  rate  constants  are  expressed  in  the  modified- 
Arrhenius  form: 

k(T)  «  AT2  exp(C/T) 

where  T  is  the  temperature  (°K)  and  A,  B  and  C  are  constants. 


(41) 


k.1.5  Enthalpy  and  Specific  Heats 

The  enthalpy  of  each  species  is  expressed  as 

-hoj  +  VT) 

where  h0j  is  the  heat  of  formation  of  species  J  at  0°K  and  hj  is  the 
temperature  dependent  (sensible)  enthalpy.  These  are  available  for  a  wide 
variety  of  species  in  the  JANAF  tables  (5l].  For  computational  convenience, 
the  temperature  dependent  enthalpy  has  been  fit  to  a  power  series  in  tempera¬ 
ture  as 

hj(T)  -hji*  v  +  v2t - Vs  •  <k3) 


The  specific  heat  at  constant  pressure  for  each  of  the  species  j  is 
calculated  from  their  enthalpies  hj  by 


dT  1  P 


(U4) 


U.2  Input  Parameters  for  the  Hydrogen-Oxygen-Nitrogen  System. 

The  chemical  kinetic  rate  scheme  (50]  used  involves  the  eight  reactive 
species  H2,  02,  0,  H,  OH,  H02,  H202,  H20  and  diluent  which  is  considered  to 
be  N2»  The  reaction  scheme  and  the  constants  in  the  reaction  rate  expression 
(E-q*  hi)  are  presented  in  Table  I.  This  scheme  has  been  extensively  tested 
against  experimental  data  (50,52]  and  shown  to  give  good  results.  Burks  and 
Oran  (50]  showed  that  the  results  computed  with  this  scheme  conpared  very 
well  with  experimentally  observed  induction  times,  second  explosion  limits 
and  the  temporal  behavior  of  reactive  species.  Oran,  et  al.  [52]  have  shown 
that  the  scheme  gives  good  results  when  coupled  with  a  fluid  dynamic  model  in 
the  simulation  of  the  conditions  behind  a  reflected  shock. 
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Ifcble  I.  H2-02  Elementary  Reactive  Mechanism 


ki  * 

AT®  exp  (-C/T)(a> 

Reaction 

A(b) 

3 

H  +  HO  t  0  +  H2 

l.i«0<-lU) 

3.00(-H*) 

1.00 

1.00 

3.50(+03) 

l*.l*8(+03) 

H  +  H02  t  H2  +  02 

l*.20(-ll) 

9«10(-ll) 

0.00 

0.00 

3.50(+02) 

2.9K+0U) 

H  +  HOj  t  HO  +  HO 

1*.20{-10) 

2«G0(-ll) 

0.00 

0.00 

9.50(+02) 

2.02(+0U) 

H  +  H02  X  0  +  H20 

8.30(-ll) 

1.75(-12) 

0.00 

0.1*5 

5.00(+02) 

2.81*(+0l<) 

h  +  h2o2  t  ho2  +  h2 

2.80(-12) 

1.20(-12) 

0.00 

0.00 

1.90(+03) 

9.1*0(+03) 

h  +  h2o2  i  ho  +  h2o 

5.28(-10) 

3.99(-10) 

0.00 

0.00 

l*.50(+03) 

l*.05(+0l*) 

HO  +  H2  t  H  +  H20 

1.83(-15) 

1.79(-lk) 

1.30 

1.20 

1.81*(+03) 

9.6l(+03) 

HO  +  HO  t  H2  +  02 

1.09(-13) 

2.82(-ll) 

0.26 

0.00 

l.l*7(+0l<) 

2.1*2(+0l*) 

HO  +  HO  t  0  +  H20 

l.OOC-16) 

3.20(-15) 

1.30 

1.16 

0.00(+00) 

8.77(+03) 

HO  +  H02  t  H20  +  02 

8.30(-ll) 

2.38(-10) 

0.00 

0.17 

5.03(+02) 

3.691+01*) 

HO  +  H20  t  H02  +  H2 

1.70(-11) 

k.70(-U) 

0.00 

0.00 

9.10(+02) 

1.65(+OU) 

HO  +  H2  X  HO  +  H20 

1. 20( -12) 
1.33(-ll*) 

0.00 

0.1*3 

9«l*l(+03) 

3.62(+0l*) 

ho2  +  ho2  %  h2o2  +  o2 

3.00(-U) 

1.57(-09) 

0.00 

-0.38 

5*00(+02) 

2.20(+0l*) 

References 


(561 

156) 

(56) 

156) 

(56| 

(561 

(57) 

kr  *  kf/!Cc 

(56) 

(56) 

(56) 

kr  -  VKc 

(58) 

(58| 

kf  -  k  /K 

(59)  C 

(581 

kr  -  kf/Kc 

(60) 

kr  =  kf/Kc 

(56) 

(56| 

(59) 

kr  *  kf/'-‘.c 

(57) 

kr  *  VKC 
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Table  I.  H?-02  Elementary  Reactive  Mechanism 
(Continued) 


Reaction 

kj.  «  AT® 

A<b) 

exp  (-C/T)(a) 

B 

c(b) 

References 

0  +  HO  X  H  +  0, 

2.72(-12) 

0.28 

-8.10(+01) 

kf  -  k  /K 

3.70(-10) 

0.00 

8.1*5(+03) 

f(561r  C 

0  +  HO,  $  HO  +  0, 

8.32(-ll) 

0.00 

5.03(+02) 

[60| 

2.20(-ll) 

0.18 

2.82(+04) 

kr  -  kf/Kc 

0  +  H?0?  X  H?0  +  0? 

1.  U0(  -12) 

0.00 

2.12(+03) 

(57) 

5.70(-lU) 

0.52 

1*.48(+0U) 

kr  -  kf/Kc 

0  +  H,0,  X  HO  +  HO, 

1.40(-12) 

0.00 

2.13(+03) 

1 57  ] 

2.07(-15) 

0.64 

8.23(+03) 

kr  “  VKc 

H  +  H  +  MtH,  +  M 

1.80(-30) 

-1.00 

0.00(+00) 

[56] 

3.70(-10) 

0.00 

4.83(-04) 

[56] 

H  +  HO  +  M  X  H,0  +  M 

6.20(-26) 

-2.00 

0.00(+00) 

[56] 

5.80(-09) 

0.00 

5.29(+04) 

[561 

H  +  0,  +  M  X  HO,  +  M 

U.lM-33) 

0.00 

-5.00(+02) 

[56] 

3.50(-09) 

0.00 

2.30(+04) 

[561 

HO  +  HO  +  M  X  H,0,  +  M 

2.50(-33) 

0.00 

-2.55(+03) 

(56] 

2.00(-07) 

0.00 

2.29(+04) 

[56] 

O+H+MtHO+M 

8.28(-29) 

-1.00 

0.00(+00) 

[6li 

2.33(-10) 

0.21 

5.10(+04) 

kr  =  kf/Kc 

0  +  HO  +  M  X  HO,  +  M 

2.80(-3l) 

0.00 

0.00(+00) 

[61] 

l.lO(-OU) 

-0.1*3 

3.22(+04) 

kr  -  VKc 

0  +  0+  M  tO,  +  M 

5*  20( —35 ) 

0.00 

-9.00(+02) 

[561 

3.00(-06) 

-1.00 

5.94(+04) 

(561 

(a)  Bimolecular  reaction  rate  constants  are  given  in  units  of  cm3/ (molecule  sec). 
Termolecular  reaction  rate  constants  are  given  in  units  of  cm6/ (molecule 2  sec). 

(b)  Exponentials  to  the  base  10  are  given  in  parenthesis;  i.e. ,  1.00(-10)  ■ 

1.00  x  10‘10. 
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For  the  calculation  of  the  thermal  conductivity,  the  binary  diffusion 
coefficients  and  the  thermal  diffusion  ratios,  we  need  the  collision 
(l  l)#  (l  2 (2  2)* 

integrals  $r  *  '  ,  Sr  *  and  Or  *  '  .  These  depend  on  the  interoolecular 
potential  function  chosen.  For  the  results  presented  in  the  next  section,  we 
have  used  the  Lennard-Jones  (6-12)  potential  function  which  adequately  de¬ 
scribes  the  interaction  of  non-polar  molecules.  This  potential  function  can 
be  written  as 

♦(r)  *  ke[( f)12  -  (r)']  (l,5) 

where  e  is  the  depth  of  the  potential  well  (the  maximum  energy  of  attraction) 
and  0  is  the  collision  diameter  for  low  energy  collisions.  The  collision 
integrals  for  this  potential  function  have  been  calculated  and  tabulated  in 
many  independent  investigations  (e.g.  Kef.  [46],  [53]).  We  have  U3ed  the 
tabulations  of  Klein  et  al  (531 •  There  is  some  uncertainty  in  the  values  of 
the  parameters  a  and  c  in  Eq.  (45).  The  parameters  we  have  used  are  pre¬ 
sented  in  Table  II  and,  except  for  <jjj,  anc*  eh20  are  tllose 

given  by  Svehla  [54].  Since  H20  is  a  polar  molecule  the  Lennard-Jones  poten¬ 
tial  function  is  a  poor  approximation.  A  Stnckmayer-type  potential  function 
may  be  more  applicable  to  polar  molecules  [54].  However  we  have  assumed 
that  the  Lennard-Jones  potential  function  is  valid  for  H20,  but  have  used  the 
values  given  by  Dixon-Levis  [33]  for  q  and  e^Q  (and  for  o^) 
which  he  has  found  to  be  satisfactory  in  extensive  investigations  of  the 
H2-O2  system. 

The  ordinary  diffusion  coefficients  have  been  calculated  using  Eq.  (33), 
where  the  values  of  A ^  and  have  been  taken  from  the  data  given 
by  Mason  and  Marrero  [49].  However,  no  such  data  is  available  for  many  pairs 
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Tkble  II.  Lennard-Jones  (12:6)  Potential  Parameters 


Species 

e^k  (°k) 

o(lO“8  cm) 

H 

37 

3.5 

0 

106.7 

3.05 

H2 

59.7 

2.827 

OH 

79.8 

3.147 

h2o 

260.0 

2.8 

°2 

106.7 

3.467 

HO  2 

106.7 

3.467 

h2o2 

289.3 

4.196 

n2 

71.4 

3.798 

of  species  occurring  in  the  hydrogen-oxygen-nitrogen  system.  For  these  pairs 

of  species  the  diffusion  coefficients  have  been  estimated  by  noting  from  Eq. 

(34)  that  if  two  pairs  of  species  (j,k)  and  (l,m)  have  similar  values  for  the 
(l  l)* 

ft  *  '  term,  then  their  binary  diffusion  coeffients  can  be  related  by 


where 


M) 


*  MJ  +  \ 


(47) 


and  cfjk  *  ■  g  --k  .  (48) 

However,  this  procedure  only  modifies  the  temperature-independent  term  A  ,. 

When  for  the  two  pairs  of  species  is  very  different  the  temperature 

dependence  (and  hence  V  will  also  have  to  be  modified.  The  values  for 
Ajk  and  Bjk  for  all  the  pairs  of  species  occurring  in  the  hydrogen-oxygen- 

nitrogen  system  are  presented  in  Table  III.  Using  Eq.  (33)  rather  than 
Eq.  (34)  avoids  using  tables  to  evaluate  the  collision  integral  ft^*1^  for 

each  pair  of  species  at  each  temperature. 

For  the  thermal  diffusion  ratio  (Eq.  38),  the  ratio  of  colli3on  inte- 

* 

grals  C..  has  been  assumed  to  be  unity.  Although  the  collision  integrals 
<3K 

(l  l)*  (l  2)* 

ft  ’  and  ft  ’  '  are  different  from  unity  and  vary  with  tenperature,  their 
ratio  is  close  to  unity  and  hardly  varies  with  temperature.  The  collision 
integrals  and  their  ratio  for  a  range  of  temperatures  (300K  -  3000K)  for  the 
H-N2  pair  is  given  in  Table  IV. 
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Note:  For  each  pair  of  species,  the  upper  term  is  A.^  and  the  lower  term  is  B.^; 
e.g.  Dh_0  =  6.3  x  1017  T°*728/N. 


Table  IV.  The  Ratio  C*fcas  a  Function  of  the  Reduced 
Tenperature  T*  for  H-N2 


T(°K) 


» 

T 


n 


J* 


nCl>2> 

Jk 


'J* 


300 

500 

1000 

2000 


5.837 

9*728 

19.456 

38.912 

58.368 


0.T6U58 

0.70387 

0.63215 

0.56789 


0.81740 

0.74570 

O.66696 

0.59886 

0.56233 


0.9354 

0.9439 

0.9478 

0.9483 

0.9480 


3000 


0.53309 


5.  IGNITION,  QUENCHING  AND  PROPAGATION  OF  FLAMES  IN  HYDROGEN-OXYGEN-NITROGEN 
MIXTURES. 

Hie  model  described  above  has  been  used  to  study  a  variety  of  problems. 
These  include  the  calibration  of  a  theoretical  model  for  the  ignition  of  pre¬ 
mixed  gases  [6,7 1 »  calculations  of  flammability  limits  (8]  and  the  estimation 
of  burning  velocities  [ 9l •  In  this  section  we  first  discuss  some  of  the 
results  obtained  in  a  study  of  the  ignition  and  quenching  of  a 
2:1:10/H2:02:N2  mixture  [7] •  This  is  an  inherently  timeniependent  problem 
and  we  can  take  advantage  of  this  property  of  the  model.  The  second  problem 
discussed  is  a  simulation  of  a  propagating  flame  in  a  stoichiometric  H2«<dr 
mixture.  The  steady  state  burning  velocity,  that  is,  the  velocity  of  the 
flame  relative  to  that  of  the  unburned  gases,  has  been  estimated  from  the 
time-dependent  simulations.  In  both  cases,  the  initial  temperatures  and 
pressures  of  the  unburned  gases  were  298K  and  1  atm,  respectively. 

5.1  Minimum  Ignition  Energies  and  Quenching  Distances 

The  model  was  configured  in  spherical  geometry  with  an  open  boundary  at 
one  end  to  simulate  a  very  large  system.  Energy  was  deposited  at  the  center 
with  a  radius  of  deposition,  RQ.  Results  from  a  typical  calculation  are 
presented  in  Fig.  1.  The  figure  depicts  the  time  history  of  the  tenquerature 
profile  after  1+  mJ  of  energy  is  deposited  linearly  over  a  period  of  0.1  ms. 
Even  after  the  energy  deposition  is  stopped,  the  central  temperature  con¬ 
tinues  to  increase  due  to  the  heat  released  in  chemical  reactions.  With 
time,  however,  the  tenqperature  near  the  center  decreases  and  the  temperature 
away  from  the  center  increases  due  to  diffusive  transport.  By  1*.5  ms,  we  see 
that  the  temperature  distribution  exhibits  a  "flame"  temperature  profile. 
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TEMPERATURE  l*K) 


RADIUS  (cm) 

Fig.  1  —Time  history  of  the  temperature  profile  in  a  H2:O2:N2/2:l:10  mixture 
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When  the  amount  of  energy  deposited,  was  reduced  to  3  mJ,  the  tempera¬ 
ture  distribution  did  not  develop  into  a  flame  temperature  profile. 

By  repeating  the  confutations  for  different  values  of  EQ,  a  bound  for 
the  minimum  ignition  energy  for  that  particular  radius  was  obtained.  Similar 
calculations  were  performed  for  different  values  of  the  radius  of  deposi¬ 
tion,  Rq.  The  results  of  such  investigations  are  shown  in  Fig  2.  A  propa¬ 
gating  flame  results  when  3.8  aJ  of  energy  is  deposited  in  a  sphere  with  a 
radius  of  0.1  cm.  However  if  the  same  amount  of  energy  is  deposited  in  a 
sphere  of  smaller  radius,  the  rate  of  heat  liberation  is  insufficient  to 
compensate  for  the  rate  of  heat  loss  and  consequently  there  is  no  ignition. 
This  radius,  0.1  cm,  is  the  "quench-radius"  for  this  particular  mixture.  For 
radii  slightly  larger  than  the  quench-radius,  the  minimum  ignition  energy  is 
almost  constant  and  for  larger  radii  (larger  than  0.11  cm)  the  minimum  igni¬ 
tion  energy  increases  rapidly  with  increasing  radii.  Therefore  for  the  sys¬ 
tem  under  study*  the  absolute  minimum  energy  is  about  3.7  mJ.  These  observa¬ 
tions  are  in  qualitative  agreement  with  those  of  Lewis  and  von  Elbe  (55)* 
Quantitative  comparisons  are  not  possible  since  the  conposition  of  the  mix¬ 
ture  and  the  time  for  energr  deposition  are  different. 

5.2  Flame  Propagation  in  Stoichiometric  H?-Air  Mixtures  at  298K  and  1  atm. 

The  model  described  above  can  also  be  used  to  study  the  propagation  of 
flames  in  premixed  gases.  As  shown  in  Fig.  1,  if  sufficient  energy  is  depos¬ 
ited  at  the  center  of  a  sphere  and  computations  are  carried  out  for  a  suffi¬ 
ciently  long  time,  the  temperature  distribution  attains  a  typical  "flame 
profile".  However,  this  method  is  an  expensive  way  to  generate  a  propagating 
flame  since  so  much  time  is  spent  in  establishing  and  overcoming  the  initial 
conditions.  Another  method  for  initializing  the  problem  quickly  is  to  start 
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the  computations  with  a  good  guess  for  the  temperature  and  species  profile  in 
a  region  behind  a  flame  front.  The  closer  the  initial  profiles  are  to  the 
final,  steady  state  flame  profile,  the  sooner  the  initial  conditions  relax  to 
the  steady  propagating  flame. 

For  the  results  discussed  below,  the  model  was  set  up  in  Cartesian  geom¬ 
etry  with  one  end  closed  (preventing  any  gas  flow  through  it)  and  the  other 
open  to  the  atmosphere.  At  the  closed  end,  a  first  approximation  to  the 
tenperature  and  major  species  profile  was  set  up  as  initial  conditions. 

In  order  to  evaluate  the  importance  of  the  thermal  diffusion  in  the 

propagation  of  flames  we  have  studied  two  cases.  For  the  first  case  the 

effects  of  thermal  diffusion  have  not  been  considered,  that  is  the  source 

terms  due  to  thermal  diffusion  in  Eq.  (8)  have  been  neglected.  For  case  2, 

T 

the  thermal  diffusion  ratio  K  j  of  all  species  have  been  calculated  using 
Eq.  (38).  However,  in  both  cases  the  heat  flux  due  to  concentration  gradi¬ 
ents  (the  Dufour  effect)  has  not  been  considered  since  it  is  expect  to  be 
negligible  even  when  thermal  diffusion  is  not  negligible  (lO) . 

Case  1:  Thermal  Diffusion  Neglected. 

The  calculations  were  initialized  as  described  above  and  soon  a  propaga¬ 
ting  flame  developed  as  can  be  seen  in  Fig.  3,  where  the  temperature  profiles 
at  15  us  and  20  us  (from  the  start  of  the  computations)  are  shown.  In  order 
to  determine  the  flame  speed,  a  criterion  for  the  location  of  the  flame  is 
required  since  the  flame  is  of  finite  thickness.  One  method  is  to  choose  an 
arbitrary  value  of  temperature  in  the  high  gradient  region  to  define  the 
location  of  the  flame.  The  movement  of  the  location  of  this  value  in  time 
and  space  gives  the  rate  of  propagation  of  the  flame.  The  location  of  the 
flame  (defined  by  900  K)  is  presented  as  a  function  of  time  in  Fig.  4.  The 
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Fig.  3  —  Temperature  profiles  in  a  propagating  flame 
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TIME  (ms) 

Fig.  4  —  Time  history  of  the  location  of  a  flame 
in  a  stoichiometric  hydiogen-air  mixture 
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!  slope  of  this  curve  gives  the  flame  velocity  for  an  inertial  observer  and  it 

is  shown  as  curve  (a)  in  Pig.  5.  The  flame  propagates  initially  at  a  veloc¬ 
ity  of  about  13  m/s  but  by  20  ys  (from  the  start  of  the  confutations)  it 
attains  a  nearly  constant  value  of  10  m/s. 

A  quantity  of  considerable  practical  interest  is  the  burning  velocity 
•  which  is  defined  as  the  rate  at  which  the  flame  consumes  the  reactants,  or 

equally  well  as  the  velocity  of  the  flame  relative  to  that  of  the  unburned 

t 

gases.  In  Fig.  6  the  fluid  velocity  profile  in  the  system  is  depicted  at  a 
particular  time.  The  velocity  is  zero  at  one  end  since  that  end  is  a  closed 
boundary.  The  flow  velocity  increases  across  the  flame  and  is  maxinum  at  the 
open  boundary.  An  estimate  for  the  burning  velocity  can  be  obtained  by  sub¬ 
tracting  the  velocity  of  the  unburned  gases  from  the  flame  velocity.  How¬ 
ever,  there  i3  3ome  uncertainty  in  the  evaluation  of  the  velocity  of  the 
unburned  gases  since  the  location  of  the  flame  front  itself  is  arbitrary. 

The  temperature  profile  across  the  system  has  also  been  depicted  in  Fig.  6. 
From  this  figure  we  see  that  an  upper  estimate  for  the  velocity  of  the  un¬ 
burned  gases  is  the  fluid  velocity  at  the  open  boundary.  The  flow  velocity 
at  the  open-boundary  (which  is  also  the  rate  of  movement  of  the  open  boundary 
because  of  the  Lagrangian  coordinate  system  used)  is  shown  in  Fig.  5  as  curve 
(b).  Corresponding  to  the  "steady"  flame  velocity  of  10  m/s,  the  velocity  of 
i  the  unburned  gp  es  is  about  7*85  m/s.  This  gives  a  burning  velocity  of  2.15 

'  1  ,  m/s  for  t't  ■  stoic  iiometrio  hydrogen-air  flame  studied. 

i  ; 

j  There  are  small  oscillations  in  the  flame  velocity  shown  in  Fig.  5(a) 

i  and  these  lead  to  uncertainties  in  the  estimate  of  the  "steady"  flame  veloc- 

I 

1  ity.  The  amplitude  of  the  oscillations  decreases  when  the  spatial  resolution 

j  is  increased.  This  can  be  seen  in  Fig.  7  where  the  time  history  of  the  flame 

j  : 
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--  Time  history  of  the  velocity  of  (a)  the  flume  and 
(b)  the  fluid  nt  the  open  boundary 
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Fig.  7  —  Effect  of  spatial  resolution  on  the  time  history  of  the  flame  velocity 
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velocity  computed  using  100  grid  points  (cells)  is  compared  to  that  using  150 
grid  points.  Since  the  computations  are  initialized  using  the  uniformly 
spaced  grid,  when  the  number  of  grid  points  is  increased  the  spatial  resolu¬ 
tion  increases  all  across  the  system  even  in  the  region  veil  ahead  of  the 
flame  front. 


Case  2.  Effect  of  Thermal  Diffusion  on  Burning  Velocity 

In  this  case  the  thermal  diffusion  source  terms  in  Eq..  (8)  were  included 
in  the  calculation  of  the  diffusion  velocities.  The  flame  velocity  and  the 
velocity  of  the  unburned  gases  were  determined  as  described  above  for  case  1. 
The  time  history  of  these  velocities  have  been  compared  to  those  of  case  1  in 
Fig.  8.  The  effect  of  thermal  diffusion  is  to  reduce  both  the  flame  velocity 
and  the  velocity  of  the  unburned  gases.  The  "steady"  flame  velocity  (corre¬ 
sponding  to  20  ys  from  the  start  of  the  calculations)  decreases  from  10  m/s 
to  about  9*65  m/s  while  the  corresponding  velocity  of  the  unburned  gases 
decreases  from  7.85  m/s  to  about  7.65  m/s.  The  net  effect  is  to  lover  the 
burning  velocity  of  the  stoichiometric  hydrogen-air  mixture  to  about  2  m/s. 
The  reasons  for  these  effects  are  examined  below. 

The  source  terms  considered  in  Eq.  (7)  for  the  calculation  of  diffusion 
velocities  are: 

(a)  due  to  concentration  gradients  (ordinary  diffusion) 

s0D  -  v(^)  m 

(b)  due  to  tenperature  gradients  (thermal  diffusion) 

S^d  -  Kj  .  (50) 

The  spatial  variations  of  these  two  terms  have  been  calculated  for  the  steady 


41 


propagating  flame  coresponding  to  20  us  in  Fig.  8.  The  temperature  profile 
across  the  system  at  20  us  is  given  in  Fig.  9*  For  convenience,  the  abscissa 
is  given  in  cell  numbers.  In  Fig.  9,  the  spatial  variation  of  the  source 
terms  SQD  and  STD  have  been  compared  for  the  species  02  and  N2.  We 
see  that  for  02,  the  term  is  opposite  in  sign  to  Sq^.  A  negative 
value  for  Srpp  implies  that  the  oxygen  molecules  are  diffusing  from  left 
to  right,  that  is,  from  the  hotter  region  to  the  cooler  region.  Since  the 
sum  of  the  two  source  terms  determines  the  diffusion  velocity,  the  effect  of 
thermal  diffusion  is  to  slow  the  diffusion  of  oxygen  molecules  into  the 
higher  temperature  region  of  the  flame.  The  effect  of  this  on  the  burning 
velocity  is  not  obvious  since  it  depends  on  the  details  of  the  reaction 
mechanism  and  the  distribution  of  other  radicals.  For  the  species  N2,  the 
term  STp  is  comparable  to  the  term  SQD  and  the  effect  of  thermal 
diffusion  is  to  enhance  the  movement  of  N2  from  the  hotter  region  of  the 
flame  to  the  cooler  region. 

In  Fig.  10,  we  compare  the  terms  SqD  for  the  species  H2  and  H.  For 
H2,  the  effect  of  thermal  diffusion  is  similar  to  that  of  ordinary  diffusion 
and  it  enhances  the  movement  of  H2  from  the  cooler  region  to  the  hotter 
region.  The  term  Sgp  for  the  species  H  is  particularly  interesting.  It 
is  observed  that  the  H  atoms  in  the  region  up  to  cell  number  75  (i.e.  T  > 
lUOO  K  have  a  tendency  to  diffuse  into  the  hotter  region  while  those  after 
cell  number  75  (T  <  1400  K)  diffuse  into  the  cooler  region.  The  diffusion  of 
H  atoms  into  the  cooler  region  is  very  important  in  the  propagation  of  the 
flame.  The  effect  of  thermal  diffusion  is  to  delay  and  decrease  the  diffu¬ 
sion  of  H  atoms  into  the  cooler  region.  Since  the  H  atoms  are  produced  pri¬ 
marily  in  the  high  tenperature  region  (in  fact  here  the  concentration  of  H 
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Fig.  9  —  Comparison  of  the  spatial  variation  of  the  source  terms  due  to  ordinary  diffusion 
(Sod)  dermal  diffusion  (S»j<d)  f°r  the  species  O2  and  N2.  The  temperature  profile 
across  the  system  has  also  been  depicted. 
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Fig.  10  —  Comparison  of  the  spatial  variation  of  the  source  terms  due  to  ordinary  diffusion 
(Sod)  an<*  thermal  diffusion  (S-pjj)  for  the  species  H  and  Hq 
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attains  its  maxima  at  cell  number  75),  thermal  diffusion  of  H  can  signifi¬ 
cantly  retard  the  rate  of  propagation  of  the  flame.  This  may  be  the  min 
reason  for  the  observed  reduction  in  the  burning  velocity  when  the  effects  of 
therml  diffusion  vere  included. 

The  term  STD  for  the  other  species  in  the  system  is  small  con^ared 
to  the  term  S0D.  In  Fig.  11  ve  also  see  that  the  source  term  due  to 
therml  diffusion  of  oxygen  atoms  is  significantly  lesr*.;  than  that  due  to 

hydrogen  atoms. 
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6.  SUMMARY  AND  CONCLUSIONS 


This  report  described  a  one-dimensional,  time-dependent ,  Iagrangian 
model  for  predicting  the  initiation,  propagation  and  quenching  of  laminar 
flames  in  pre-mixed  gases.  A  number  of  new  approaches  and  algorithms  used  in 
the  model  have  been  discussed.  An  asymptotic  coupling  method  was  used  in 
conjunction  with  timestep  splitting  to  couple  the  various  physical  and  chemi¬ 
cal  processes.  This  approach  allowed  the  use  of  entirely  different  algo¬ 
rithms  for  the  processes  represented  by  different  terms  in  the  conservation 
equations.  This  made  the  computing  procedure  very  efficient  and  inexpen¬ 
sive. 

The  result  from  several  calculations  using  the  model  have  been  pre¬ 
sented.  The  first  is  a  flame  initiation  and  minimum  ignition  energy  study  of 
a  hydrogen,  oxygen,  nitrogen  mixture.  This  used  the  time-dependent  property 
of  the  model  and  the  oberv&tions  were  in  qualitative  agreement  with  experi¬ 
mental  data.  Quantitative  comparisons  were  not  possible  since  the  composi¬ 
tion  of  the  mixture  and  the  time  for  energy  deposition  were  different.  The 
second  problem  discussed  was  an  estimation  of  the  steady  state  burning  veloc¬ 
ity  of  a  stoichiometric  hydrogen-air  mixture.  The  burning  velocity  predicted 
by  the  model  was  on  the  lower  end  of  the  observed  experimental  range.  We 
note,  however,  that  there  are  several  factors  which  would  tend  to  increase 
this  estimated  value.  First,  the  calculated  burning  velocity  would  be  larger 
if  the  velocity  of  the  unburned  gas  used  in  this  calculation  was  closer  to 
the  flame  front  instead  of  at  the  open  boundary.  It  is  not  clear  where  it 
should  be  evaluated  in  a  time-dependent  calculation.  Second,  the  uncertain¬ 
ties  in  the  values  of  some  of  the  diffusion  coefficients  are  about  30^.  The 


burning  velocity  was  found  to  be  very  sensitive  to  certain  diffusion 
coefficients#  Because  of  the  uncertainty  in  the  experimental  values  of  the 
diffusion  coefficients,  there  may  not  be  any  advantage  in  using  the  empirical 
formula  over  the  expression  from  kinetic  theory. 

The  model  also  predicts  that  the  effect  of  thermal  diffusion  is  to  lower 
the  burning  velocity  by  about  1%,  This  observation  is  in  agreement  with  that 
of  Dixon-Lewis  [32]  who  found  that  neglecting  thermal  diffusion  of  hydrogen 
atoms  increases  the  burning  velocity  by  5  or  6%.  For  the  hydrogen  atom,  the 
thermal  diffusion  source  term  was  found  to  be  small  compared  to  the  source 
term  due  to  ordinary  (concentration)  diffusion.  But  because  of  the  complex 
chemistry  occurring  in  flames,  a  small  difference  in  the  hydrogen  source  term 
causes  a  significant  difference  in  the  burning  velocity.  This  again  empha¬ 
sizes  the  need  for  more  accurate  estimates  of  the  diffusion  coefficients. 

In  conclusion,  the  model  can  be  used  confidently  to  predict  the  igni¬ 
tion,  propagation  and  quenching  of  laminar  flames  in  pre-mixed  gases  provided 
there  is  satisfactory  data  on  the  chemical  kinetics  and  the  diffusive  trans¬ 
port  coefficients. 
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